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ABSTRACT 

Asymptotic-induced methods are presented for the numerical solution of hyperbolic con- 
servation laws with or without viscosity. The methods consist of multiple stages. The first 
stage is to obtain a first approximation by using a first-order method, such as the Godunov 
scheme. Subsequent stages of the method involve solving internal-layer problems identified 
by using techniques derived via asymptotics. Finally, a residual correction increases the 
accuracy of the scheme. The method is derived and justified with singular perturbation 
techniques. 
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1 Introduction 


The combination of asymptotic and numerical analyses provides improved accuracy for 
multiple-scales problems 2 . Many physical problems have multiple scales; a typical situa- 
tion occurs when physics on the fastest scale induces narrow regions where the variation in 
the solution is large. Such regions are called boundary layers or transition layers , depending 
on whether they are near a boundary or inside the interior of the domain. Examples of 
such situations are laminar flow of a slightly viscous fluid or combustion with high activa- 
tion energy. Classical schemes applied to these types of situation generally fail to correctly 
describe the behavior inside the layers. This difficulty is overcome by developing numerical 
methods based on the asymptotic analysis of the multiple-scales problems. We demonstrate 
how asymptotic analysis and numerical analysis can interact to achieve the goal of a highly 
accurate solution method. 

The numerical techniques developed in this paper are designed to solve two problems: 
hyperbolic conservation problems with or without viscosity. Also, we treat three different 
types of singularity arising in inviscid and viscous conservation laws, using domain decom- 
position: specifically, we solve problems that have shocks , weak singularities t and interaction 
of singularities. We use an asymptotic analysis to identify the appropriate scalings as well 
as modified problems in the transition layers. The domain decomposition of our numerical 
method is then based on the criterion given by the asymptotic analysis. The numerical algo- 
rithms involve several stages to solve each modified problem in the layers. Finally, we treat 
weak singularities with a residual correction. The residual correction also identifies incorrect 
asymptotic preconditioning and controls the numerical accuracy of the solution. 

A number of recent papers present new numerical methods that are induced by an asymp- 
totic analysis. For example, basis functions motivated by asymptotics of boundary layers 
are developed in [3]. Many of the basic ideas relating to asymptotic analysis and numerical 
methods that use domain decomposition are found in [2]. These ideas were incorporated into 
a parallel numerical method in [15]. Specific applications to problems involving conservation 
laws have been developed in [1]. 

On the other hand, some recent papers used the singular perturbation approach to in- 
vestigate the effect of viscosity on conservation laws. Sharp estimates can be found, for 
example, in [10] and its references. The dynamics of internal layers are studied in [6]. Also, 
uniform approximation based on the matching asymptotic technique are given in [8] and [7]. 

This paper is a further stage toward the synthesis of asymptotic analysis and numerical 
methods applied to systems of viscous conservation laws (seen as a singular perturbation 
problem). We begin with an asymptotic analysis of the problem in Section 2. This anal- 
ysis is expanded and combined with numerical techniques in Section 3. We end with a 
demonstration of the techniques on computational fluid dynamics problems in Section 4. 

2 Asymptotic Analysis 

In this section we present the problem and briefly review the asymptotic analysis that is used 
in the method. The result is a uniform approximation based on the asymptotic technique of 

2 In this manuscript we shall solve all important problems of current interest, except possibly inflation. 
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matching [5]. 

Asymptotic analysis that is strongly coupled with the numerical analysis is delayed until 
later in the paper. 


2.1 Setting of the Problem and Regular Expansion 

Consider the Cauchy problem 

I f + T,nV) = (Wf ) for (x,i)sn 

I U(x, 0) — V(x ) for igE. 


( 1 ) 


Here the solution U : fi — * IR n is a vector-valued function , the domain is 0 = IRx]0,T[, 
and e << 1 is a small parameter. 

We assume that V is piecewise smooth. We also assume that F and P are smooth 
functions of U . We suppose that P is a suitable viscosity matrix [4] for the shocks of the 
following associated inviscid problem: 


r f + £*W = o for (x,t)en 
k U{x, 0) = V(x) for x £ IR. 


( 2 ) 


That is, a shock wave solution to (2) can be obtained as a limit of progressive wave solutions 
of (1). Problem (1) is a parabolic-hyperbolic singular perturbation problem driven by (2). 
One easily obtains the regular expansion 

pouter = U 0 + eU l + e 2^2 + (3) 


that is a priori valid outside the neighborhood of the singularities of the solution to (2). We 
substitute U°" teT in the differential equation of (1) and use identification in e to obtain that 
U° must be a solution of (2). We also find that U 1 must be a solution of the following linear 
hyperbolic problem: 

TT + fi (DF^W 1 ) = £ (P(U 0 )^) for (.,!)£» 

C/ 1 (x,0) = 0 for x £ IR, 

where DF denotes the Jacobian matrix of F. In general, the equation for the term U l with 
coefficient e* is the solution to the analogous differential equation 

^ + A ( D F(f/»)^) = RHS(U’,j < i), 

with a right-hand side that depends on the previous terms. The inviscid problem (2) has 
many weak solutions; we shall uniquely define U° as the analysis progresses. 

It is most common to complete U™ tei with a £/j£ ner for each singular region to obtain a 
uniformly valid approximation to a viscous problem. However, we shall develop in detail only 
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the asymptotic analysis that is required for our numerical methods. Thus we restrict our 
asymptotic analysis to complete U™ ieT with U^ CT only when the singularity is a shock layer. 
The solution is completed for weak singularities with the residual correction techniques of 
Section 3.4. This completion also results in a uniformly valid approximation. Furthermore, 
the approximation is not difficult to obtain numerically. It is easy to see that interaction of 
singularities requires a stretching in both the space and time variables: this leads to some 
parabolic layers. We shall use a numerical treatment of the interaction of singularities that 
does not require a complete analysis of these parabolic layers. 


2.2 Transition Layer That Corresponds to Shock Waves 

We assume that solutions U° of (2) are smooth except on piecewise regular curves Sk{t). It 
is assumed for t £ [io> ^i] that the S* are isolated from each other. Without loss of generality 
we may drop the subscript k, S is assumed to be smooth for t £ [io)^i]* I n particular, we 
assume no focusing of characteristics. To make this more precise, we suppose that there is 
an interval of time [to,ti] such that for each t £ [tojh] an d for each S the following limits 
exist * 

u?= lim V°,= %/(*,!), 

Let us define the change of variable: x = x — S'(t ) t and r = t. We denote U(x,t ) = 
U(x,t). We further assume that 


r)U° 8U° 

Url = lim U°r,r= E® J-, 

' i -. 0 - OT ' S — 0 + OT 


frU j d i U j 

U)i = lim ~£rz~ r i U ] r = lim -^7-, for t ,j > 0, 

*-o- dx l h 4—0+ ox x 

v°, = 0( 1), 0;, - 0(1). 

The shock layer profile will not have rapid variation, so it is appropriate to scale and 
translate only the spatial variable (and not the temporal variable). Such a transformation 
is defined by 

x — S(t) 


t = 


and r = t, 


(5) 


where we denote ?7 (£,t) = Under this transformation the differential equation of 

problem (1) becomes 


dU 

dr 


U e -^( F{ U)-S(r)u)=rf ( (p(U) d £j 


( 6 ) 


This suggests an inner expansion of the form 

pinner _ (j0 + £ (jl + £ 2^2 + (7) 

where all of the terms U i are functions of ( and r. Using this expansion in (6) and imposing 
the matching relations [5, 8] with the outer expansion C/£, uter , we derive the set of ODE 
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problems for the U'. The equation for the first term is 

+ h [F(P°) - S‘(t)U ° ) = 0, 

U° — ► Uf as £ — ► — oo 
U° -> U? as ( -* +oo 


and for the second term we have 

' ~$r ( P0° ) -V')+Tt (DF(V°) ■ U' - S\t)U') = -f?, 

- !|£/> - (U},( + . 0 as ( ^ -<x> (9) 

. II# 1 - (Ui r ( + C/° r ) 1 1 0 -t 0 as ( -< +oo, 
where ||I7|| 0 = max; =:1 .. n (|u,|) for U = )«=!■•„• More generally, we obtain 

' ~W (Ptf 0 ) ■ #‘) + f( ( DF(0° ) ■ U‘ - S‘(t)U<) = RHSiUi.j < i), 

■ ll^ 0 -s;.o ! ^iiqj'£ i -’'iu^o »{--» (10) 

. II# 0 - £5=0 - ° as f-.+oo 


where the right-hand side is a nonlinear function depending on the for j < i. Identification 
of the coefficients of the appropriate power of e in the differential equation of problem (6) 
determines RHSi. 

The temporal variable r can be considered as a parameter in the transition layer because 
the above problems require the solution of only ODEs. Also, the problems for fj x are linear 
for i > 0. Being able to treat r as a parameter is not surprising, since U° t = 0(1) and 

£^ = 0(i). 


A 

Next we discuss the existence of Uo and the uniqueness of Uo as well as the uniqueness 
of the curve S. Using a matching relation on the first spatial derivates on the terms in the 
expansion I7 ( JJJ ner , one looks for U° such that 


dU° 

dt 


0 as £ 


We integrate (8) from — oo to £ to obtain 


±oo. 


P(U ° ) 


dU° 

d( 


= H(U°) - Hffl), 


( 11 ) 


where H(U) = F(U ) — S'(t)U. Thus, the existence of a solution to (8) implies the Rankine 
and Hugoniot condition (RH) 

mn = m°). 

This means that U° and U° are critical points of the dynamical system (11). 

In addition to the RH condition, we must be able to construct the layer; thus, we assume 
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HI. 


There exists a unique trajectory for the 
dynamical system (11) from U° to U°. 


It is interesting to compare HI with the classical geometric entropy condition (GEC). For 
scalar conservation laws, Hi is equivalent to the GEC [11] [14]; however, this is not in general 
true. A number of very interesting studies have considered which conditions one can expect 
that the GEC implies HI. In particular, the GEC implies HI when P is the identity and 
the problem consists of a 2 x 2 system with the assumption of genuine nonlinearity [4, 12]. 
The existence of a viscous profile (i.e., HI) is not of the same nature as the GEC. However, 
we restrict our problems to cases where the RH condition and HI are enough to uniquely 
define U° and S (and consequently U™ teT ). 

Now we have determined U° up to a translation in the spatial variable £. In Section 2.2 
this translation will be imposed as a solvability condition on the solutions to problem (9). 


2.3 Construction of the Inner Expansion in the Shock Layer 

Here we briefly outline the construction of the solution in the shock layer. This formal theory 
for systems of conservation laws is a generalization of the construction in the case of scalar 
conservation laws [8]. For more details of this generalization, see [9]. 

Let us concentrate on the construction of the solution U 1 of problem (9). The procedure 
is first to obtain two functions and £/*(£,£) as solutions to related problems, then 

combine them such that the resulting function is a smooth solution of (9). The function U{ 
satisfies the equation from problem (9) with the left boundary condition 

U 1 - (Uq [£ + £7°,) — ► 0 as * - — oo, 


whereas U , f satisfies the equation from problem (9) with the right boundary condition 

U 1 - (£/- o y + Ul r ) - 0 as *-+oo. 

We use the solutions to these problems that are constructed using an element dU°/d£, of 
the kernel of (9). That is, we construct U l as the function 


[ + for £<0 


. U r l ((,t) + Br(t)xS§l for f >0, 

where U X V denotes the vector (ujUj)j- 1-n . Bi and B r axe some smooth vector functions 
from fi into IR that satisfy the continuity condition 

afro afro 

+ £,(i) x ^-(0,t) = V r '(0,t) + B r (t) x -gf(0,t). 

Notice that for a fixed value of t y this relation determines the vector functions up to a single 
constant. 
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One can prove [9] that U l is a smooth function iff U° satisfies the relation 

I ui ^° {( ’ t] ~ u ° w } + r ° d( ) 

r du°i 

= - DF(17 0 ) • U 1 + P(U°)-fo . 

where [■] denotes the jump across the shock. This is simply an area relation that determines 
the shift in £ for U°. Let us assume for simplicity that U and its space derivatives vanish at 
infinity. One can show that this relation is a consequence of the conservation relation 


L 


-CO dt ’ 


satisfied by solutions to (1) and (2). This relation is also be satisfied by our uniform approx- 
imation of the solution to (1), where 

U„ = (1 - H{i u ))UT T + + eU 1 ) + 0(e 2 ). 

Here ( v = (x - S(t))/e v is the intermediate variable and H is a smooth cutoff function 

f 0 if |y| > 2 

H(y) = 

{ 1 if |y| < 1. 


More precisely, we have 



0(e), 


^Ve currently have defined t/ 1 up to a kernel function of (9). Since we use only the 
first term U° in the current implementation of this method, we do not discuss how to 
determine this constant here. This construction can be pursued at any order. For example, 
the solvability condition for problem (10) with i = 2 will select the only admissible U 
solution of problem (9) [8]. 


3 Asymptotic-Induced Numerical Scheme 

Now we consider the numerical treatment of problems (1) and (2) by using the analysis in 
the previous sections. 


3.1 Localization of the Singularities 

Let {^t 1,0,1,... be the spatial discretization of step size Ax, and let {tk}k-i,... be the 
temporal discretization of step size At. For convenience we restrict the discussion on local- 
ization of singularities to scalar conservation laws. The scalar terms are distinguished from 
the system by using the notation u and / as the solution and flux function in place of U and 
F, respectively. We suppose that u is known at time f*. Starting from u(-,tk), we apply a 
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Table 1: Asymptotic Order of Residual 


Type of Zone 

Order of Residual 
ut + f{u) x 

Local Coo 
£ 

rdinates 

T 

Regular zone 

0(e) 

X 

t 

Shock layer 

0(6-) 

(x - S(t))/e 

t 

Weak singularity 

0(e 1 ' 2 ) 

(x - S(())/e‘" 

t 

Shock interaction 
with other singularities 

0(e~') 

(x — S 0 — Sit)/e 

(t - t 0 )/e 

Discontinuous with 
/ locally linear 

0(1) 

(x - So - Sit)/e^ 2 

t-t 0 

Formation of shock 

0(e~ 1/4 ) 

(s - So - 

( t - 


discretization that gives an approximation u K at the discrete spatial points {x,} at time tx, 
where K > k. Given this approximation, we discuss the problem of localization of different 
types of singularities of u at time tx- 

Localizing the singularities will use the residual obtained from the left-hand side of the 
PDE in (1). To this end, we recall in Table 1 some results of the asymptotic analysis of the 
Cauchy problem (1); see [8] and its references. Here So and Si are the constants from the 
local linearization S(t) = Sq + S\t + 0((t — t 0 ) 2 ) in the neighborhood of the singularity. 

The localization is presented for the case when u K is obtained by using the first-order 
Godunov scheme; however, this analysis requires only minor modification for several first- 
order discretizations of scalar conservation laws, provided the discretizations are conservative. 
The solution obtained via the Godunov scheme is an approximation to within 0( Ax 2 ) of the 
following Cauchy problem: 

It + £/(“) = Ax £ (G(/,u,A*,Az)f;) < 6 | h,t K ] 

■ ( 12 ) 

k u(x,t k ) given. 

We fix the ratio At/ Ax so that Ax may be treated as the parameter e in the table. 

In the more general case of a system with some first-order discretization, we obtain Wq 
as the numerical approximation of the solution to the hyperbolic problem. Let us assume 
that the function Wg is the solution to 

' f + &F(£0 = A*£(G«(f,tf,Az,At)f) (€[(*,«*-] 

< 

k U(x,t k ) given, 

where Gh is for the numerical viscosity of the specific hyperbolic scheme with Gh = 0(1). 
The detection for the numerical method is based on obtaining an approximation to the 
residual dU/dt + dF(U)/d x. 

The residual is of magnitude 0(Ax -1 ) in either a shock layer or in a zone where a shock 
interacts with some other singularities. In addition, an approximation of the viscous term 


7 



J" ?U(-,iK ) localizes some of the singularities. For example, this viscous term will be of order 
0(Ax -1 ) in a shock layer or in a zone of interaction. 

Numerically it may be difficult to determine when the residual (or viscosity) is of order 
A® -1 . This difficulty may usually be overcome by computing the numerical solution on 
two grids with different spatial spacings A®i and A® 2 . Let Ri and he approximations 
to the residuals of the solutions on these meshes. When the ratio R 1 /R 2 of the residuals is 
within some tolerance of A® 2 /A® 1 , then we treat the region as being inside a shock (or shock 
layer). It should be mentioned that asymptotic analysis and numerical experiments show 
under general conditions that strong singularities are not well localized by any approximation 
of the gradient of the solution. 

Shocks (and interactions with shocks) are easily detected and are treated in the next 
section. In Section 3.4, we shall show how some weak singularities that are much more 
difficult to localize can be absorbed in a residual correction process. 


3.2 Numerical Treatment of Shocks 


Let us first recall the numerical treatment of a simple shock for a scalar conservation law [1], 
Let ft 0 = [a, &] he the zone that includes the shock (or the shock layer). First we construct 
an approximation w to the solution of (2); next, to compute an approximation to the shock 
location S, we treat the solution near the shock (or inside the shock layer). 

The region fio is chosen so that u is an a priori valid approximation to u except possibly 
for some portion of fl 0 ; thus, w = u outside flo- The approximation w is constructed inside 
fl 0 with an extrapolation of u. This extrapolation is based on values of u outside 0 Q . We 
construct w as the connection of the extrapolated values. For t = tx fixed, we impose the 
area relation 




w risht )dx = 0 


on the discrete problem. This is presented visually in Fig. 1, where we use second-order 
extrapolation and have placed S so that the area covered by the dark grey is equal to the 
area covered by the light grey. As we discuss below, the left and right values of w are used 
for the boundary conditions for the first-order approximation of the solution with a viscous 
profile. 

Notice that the accuracy of the approximation w as a solution to the inviscid problem 
is limited only by the particular numerical techniques. Specifically, the accuracy is limited 
by the accuracy to which u is computed outside fl 0 and by the order of accuracy of the 
extrapolation. There is no limitation imposed by the asymptotic analysis on this aspect of 
the method. 

The treatment of the solution in a neighborhood of the shock (or inside the shock layer) 
depends on whether the goal is the solution of the viscous problem (1) or the inviscid problem 
(2). If the solution to inviscid problem is the goal, then we wish to minimize the viscosity. 
Viscous problems require the computation of the profile of the shock. The asymptotic 
analysis of the shock layer for a system of conservation laws is similar to the case of a 
scalar conservation law (cf. Section 2). 

For the numeric computations, the primary difference between the treatment for a scalar 
conservation law and a system is that we must verify that certain conditions are satisfied by 
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Figure 1: Approximation of shock location 


the system. Specifically, we a posteriori (during the numerical computations) verify that 

• W° and W° axe on a Rankine-Hugoniot set up to a tolerance, 

• there is a trajectory from W° to W°, and 




°Kr 

dr 


= o(i). 


The first assumption is satisfied automatically by solutions to scalar conservation laws. Now 
we shall describe in more detail the numerical treatment of the shock layer depending on 
whether one solves an inviscid or a viscous problem. 


3.2.1 Minimum Viscosity Method 

Suppose that for a particular tk , the shock location S is in [x;, Xj_|_i[, where [x,-, Xj+i[C Oo- 1° 
the most general case, S will not be a grid point, and we modify W to maintain conservation. 
This step results in modifying either the value of W-* or the value of Wf +1 (see [1] for details). 
The modification introduces the minimum of viscosity needed to use the Godunov method 
throughout the domain. This modification of Wf or is the final stage of the shock layer 
treatment when the goal is the solution to the inviscid problem. 

3.2.2 Viscous Problems 

To solve problem (1), we implement the inner layer computation as in Section 2.2. The 
matching relations used to obtain a uniform approximation are based on the numerical 
approximation of the outer expansion. Thus, we solve (8) except that W® and \V® are used 
in place of U° and U°, respectively. 
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To obtain the higher-order corrections U* for i > 0, we would solve the problems (9) and 
(10) with modified boundary conditions; however, an easy approximation can be obtained 
by computing a poor approximation Up to U 1 such that 

f W — U° for x < Xj 
* = . 

( W — U° for x > Xj. 

The numerical treatment of this problem involves using an ODE solver to obtain the 
solution to the system of ODEs (8). Notice that except when the details of the viscous 
profiles are desired, we need only obtain the values of the inner expansion at the coarse grid 
points on which we know W. Even when a higher-order numerical method is used, there are 
a large number of steps in the integration of the ODE system for the inner region. The cost 
of this integration may be reduced by using an explicit formula of approximation given by 
the asymptotic analysis in the neighborhood of the critical points (see [1]). 


3.3 Interaction of Singularities 

In regions where solutions to problem (1) contain shocks that interact with other singularities, 
we use a brute-force approach. The local scaling in both the spatial and the temporal 
variables is 

. x So t t 0 

( = — - — , r = . 

e e 

Under this transformation the PDE that governs the solution becomes 





(13) 


where U(£,t) = U(x,t). This is the equation that is solved in the regions with interactions. 
There is no expansion of the solution used here. The transformation a priori resolves all of 
the physics. This situation is reflected by all of the terms in (13) being order unity. 

The initial condition for this problem is derived by imposing C° continuity. This was 
sufficient for the problems considered in this paper. However, the computation of the residual 
is a check of the numerical accuracy of the composite scheme. A residual correction similar 
to the one in the next section could be applied to improve the result. Interface boundary 
conditions is a topic of further study. 


3.4 Residual Correction 

Second-order accuracy in both space and time is achieved through a first-order numerical 
approximation followed by a residual correction. This accuracy is obtained in regions where 
the regular expansion (3) is valid. We also discuss the residual correction in the presence 
of weak singularities, where the accuracy is increased (but not as much as for the smooth 
regions). This discussion begins with the asymptotic analysis, is followed by the choice of 
the particular discretization, and ends with the treatment of weak singularities. 
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3.4.1 Asymptotic Analysis 

The asymptotic analysis is based on the regular expansion presented in Section 2.1. Thus, 
this analysis is restricted to regions where the solution is smooth and the regular expansion 
is valid. For convenience, our analysis assumes e and At are both 0( Ax). The dependency 
of e on Ax is related to the order of the regular expansion and the hyperbolic scheme used. 

The Godunov method applied to (2) gives us W G as a first approximation to the solution 
to the parabolic problem (1); thus U = W G + 0( At, Ax). A single correction W c will be 
computed so that 

U = W G + hW c + h 2 W 2 + 0{h 3 ), (14) 

where h is a small parameter of magnitude O(At) = 0( Ax). The equation for W c is deter- 
mined with an analysis similar to that of the regular expansion. 

The viscous terms now arise from both the viscosity in (1) and from the numerical error 
in W G , resulting in the linear problem 

f ^ + l(DF{U)W c ) = R{W a ) 

( 15 ) 

{ W c (x,0) = 0 

for W c . Here, 

kR { Wa) = (PP^) - + £W) ■ < 16 > 

where v = e or v = 0, depending on whether we are solving the viscous problem or the 
inviscid problem. Notice that the terms of higher order would be included in the next term 
of the regular expansion of U . The right-hand side in the equation of (16) is the residual 
from using W G as an approximation to the solution of (1). To obtain a second-order method, 
we start with the first-order approximation W G , then compute the right-hand side of (16) to 
second-order accuracy, and finish with a first-order scheme for (15). 

The asymptotic techniques are similar in the more general case when e is not of size 
O(Ax), and/or when a different order method is used to compute W G . Consider the case 
when the numerical method is 0(Ax p , At 9 ). In this case we choose the parameter h in the 
expansion (14) to have magnitude equal to the largest of Ax p , At 9 , and v. The problems 
for the terms in the regular expansion are obtained by identification in h after substitution 
of the expansion into 

f + l F w = 4 (w|) + 

where v = 0 or e, and Pi and P 2 are terms of 0(1). We do not require specific forms for Pi and 
P 2 for the numerical method since we may use a residual in place of Ax p Pi(U) + At 9 P 2 (t/). 

3.4.2 Discretization of Correction 

Here, we choose discretizations for the numerical approximation of the residual and the 
correction scheme. Notice that these discretizations as well as the choice of the Godunov 
method for the first approximation are only examples of possible treatments. 
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Centered differences are used throughout this section; hence, the resulting finite-difference 
scheme will be a conservative discretization. 

Residual. Here we use the discretization based on centered differences and averaging. 
The temporal derivative is approximated by dWa(xi,tk+i/ 2 )/dt « 8 t W q^ 1 ^ 2 , where we use 
the notation 6 to be the centered difference operator 

f 9(V + At?) - g(r) - Arj ) 

S, - 9(V) = 2Aj; ' 

The formula for the complete residual is 

j {- [wSJ* 7 ’ + i (6,.F? +I + «*#)] (17) 

Here we use the notation Ff = F(WQ i ) and the average P^/ 2 = (P(W^ )t ) + P(Wc i+1 ))/2. 
Taylor series analysis can be used to verify that this is an 0( Ax 2 , At 2 ) approximation to the 
residual at the point ( x,t ) = (a:,-, ^+ 1 / 2 ). 

Correction Scheme. The implicit scheme based on a centered spatial derivative and a 
centered temporal derivative is used for the discretization of W c to obtain 

S t w£ 1/2 + 6 2x (DF(W£?) ■ W# 1 ) = ^ +1/2 - (18) 

This is a first-order accurate method. 

The accuracy of the discretization for the residual is balanced with the accuracy of the dis- 
cretization for the correction term. We compute Wq to first order. The residual is computed 
to second order; however, it is scaled by 1 jh for the computation of the correction. This 
results in a first-order approximation to the source term for the correction. The correction 
discretization is first order; hence the overall method is second order. 


3.4.3 Weak Singularities 


A modification of the above analysis can be used to develop the residual correction in the 
presence of weak singularities because the residual is still small in the neighborhood of a 
weak singularity (cf. Table 1). The numerical algorithm is the same for this case; however, 
the analysis (and resulting accuracy) are different. In place of (14) we have the regular 
asymptotic expansion 

U = W G + Wi + . . . , 

where Wg is the approximation of U obtained from the numerical method. The function 
Wg is an approximation to the solution of 


dU d . d 
~dt + ~di F ^ ~ V !h 



where v = e or v = 0. However, Wg may be treated as the exact solution of the modified 
equation 


dW G 

dt 


+ ^P(^g) = h RHS, 


(19) 
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where h = o(l) is defined such that RHS = 0,(1) is true a priori. So we rescale W x as: 
W x = hW c . As in Section 3.4.1, we obtain the equation for W c as 


dW c 

dt 


±(DF(Wo)-W.)=-RHS + ll 



c?Wg\ 
dx ) 


This is the equation for the correction term in the presence of weak singularities. 

The numerical discretizations used for the case of the weak singularity can be the same as 
those presented in the last section. For example, when the singularity is of order h = 0(y/e), 
then the method with the same discretizations will have accuracy 0(h 2 ) = 0(e). Both forms 
of the analysis in this section are valid with the numerical discretizations of the previous 
two sections. The numerical method after the residual correction will be more accurate 
than using the first approximation Wg alone. In addition, it is feasible to use a higher-order 
method for the discretizations of the residual and a higher-order expansion for the correction 
in regions of weak singularities. 


3.5 Numerical Algorithm 

The numerical and asymptotic analysis are combined and summarized in Algorithm 1 below. 
The outer region is discretized by using the coarsest grid. Numerical values of Wg and W c 
are determined at coarse grid points. Wg is computed over the whole grid; however, the 
coarse grid values located inside refined regions are used only in the conservation relation. 
The refined grid values are injected into these coarse grid values at the end of the time step. 
Linear interpolation is used to obtain values of the coarse grid solution between points when 
necessary. 

The numerical results presented here involve only a two-level refinement. This is because 
the current implementation is designed for problems in which shocks (or shock layers) and 
their interactions are the only strong singularities. We recall that weak singularities are 
treated by the residual correction scheme. 

The numerical solution in the shock layer is computed in the transformed coordinates 
(5) when the shock profile is not rapidly varying. When interactions exist, the coordinates 
use the same £ with the temporal scaling t = i/e. Although A£ is the same magnitude as 
Ax, local spatial scaling is used. This means that a region of size Ax in the coarse grid 
corresponds to many points in the £ -grid. For this reason, we refer to the £-grid as the 
refined grid. 

The method is adaptive. The refined grid is allowed to change shape and location as the 
solution evolves. The refined region contains an overlap region where both the inner and 
outer solutions are valid. 


4 Experiments 

In this section we demonstrate the techniques developed within this paper. The experiment in 
Section 4.1 demonstrates the residual correction. The experiments in Section 4.2 demonstrate 
the method on the isentropic gas dynamic equations. 


13 



For k = 1 


I. March from i* to tk+i on two coarse meshes with spatial discretizations Axi ^ 

A®2- 

II. Detection. 

A. On each coarse mesh compute the residual used in Table 1. (Alternatively, 
compute the viscosity.) 

B. Mark regions that should be refined. 

III. Compute Wc with the residuals of Step II.A. except in marked regions. 

A. Modify the shape of the refined region. 

IV. March the marked regions from t* to t* +1 according to the type of singularity. 

A. For a simple shock of Section 2.3, compute the solution to the ODEs as 
outlined in that section. 

B. For an interaction of singularities of Section 3.3: 

1. Form the initial condition in newly refined regions. 

2. Determine boundary conditions from linear interpolation. 

3. March inner solution Ai/(eAr) time steps. 

V. Inject the values from the internal layer to obtain the composite uniform approx- 
imation. 


Algorithm 1 Numerical algorithm 
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4.1 Solving Burgers’ Equation 

The residual correction technique is demonstrated by solving the inviscid Burgers’ equation 


du du 

w + “fc = 0 


on the domain [-0.2, 1] with a wedge initial condition. The first approximation is obtained 
with a first-order Godunov method using discretization parameters Ax = .025 and At = 
.0125. The solution at time t - . 5 is presented in Fig. 2. The increase in accuracy of 



X 

Figure 2: Exact solution to Burgers’ equation 


the residual correction over using just the Godunov is shown in Figs. 3-4. The minimum 
viscosity method is used for the numerical solution in a neighborhood of the shock. A cutoff 
function is used for the residual correction. This causes the correction to be applied only 
outside a neighborhood of the shock. The mass that might be lost or gained in this procedure 
is included in the computations for the minimum viscosity method. The large jump at the 
shock is a reflection of the error between the L 2 projection of the numerical solution with 
the exact solution. 


4.2 Solving the Isentropic Gas Dynamic Equations 

In this section the method is demonstrated on the system 

du dv 
dt dx 

dv d / 1 \ d ( dtt \ 

dt dx \u 7 / dx \dxj 

Here u is the inverse of the density and v is the velocity. These equations are obtained 
from the conservation of mass and momentum in Lagrangian coordinates assuming that u 


15 



0.2 



Figure 3: Error in whole domain with (dashed) and without (solid) residual cor- 
rection 



Figure 4: Error in corner-layer region with (dashed) and without (solid) residual 
correction 
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is equal to the pressure raised to the — l/7th power (the perfect gas law) with 7 = 2.2. The 
numerical solution in the outer region is computed with a first-order Godunov method with 
CFL n umb er At/ Ax = .2. The discretization on the scaled coordinates inside the shock layer 
is based on A£ = .1. We use a third-order Runge Kutta method to compute the viscous 
profile for a shock layer. To compute a parabolic layer, we use an explicit first-order ENO 
scheme (similar to the Godunov scheme) with the CFL condition At/A£ < .2 and stability 
condition At/A ( 2 < .25. These values are within the limits imposed for the stability of the 
finite difference method. 


4.2.1 Simple Shock 

The first experiment is a simple shock. In Fig. 5 we show the computed solution using 
various methods. The initial condition is 


u(x, 0) 
v(x, 0) 



where 

U R = 2.50, U L = .800, V R = .600, 

and Vl = Vr + 

The value for Vl is chosen using the Rankine-Hugoniot condition 

Vl - Vb 1 m - uui 

Ur -U l Vr- V t, 

The computations are run with Ax = .02. The Godunov and minimum viscosity methods 
approximate the solution to the hyperbolic problem, whereas the ODE layer (solving problem 
(8) inside the shock layer) and the parabolic-layer methods are used to compute the numerical 
approximation to the viscous problem with e = .01. 

The minimum viscosity method introduces significantly less viscosity than the Godunov 
method. As expected, the parabolic-layer and the ODE-layer methods produce numerical 
solutions that are very close to each other. 



4.2.2 Shock-Rarefaction Interaction 

The method was also used to solve a Riemann problem with a right-traveling shock and 
a left-propagating rarefaction. Thus, for the small time, we must solve an initial layer 
with a shock- rarefaction interaction. We can compute the exact solution of the inviscid 
problem. First, an analytic self-similar solution (a rarefaction emanating from the origin) to 
the inviscid isentropic gas dynamic equations is given by 

u(l , ()=7l /(, + »(£)- J/h+1 > ( 22 ) 
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Figure 5: Numerical solutions to the isentropic system, t = .4 


v(x , t ) 


2'y 1 /(l'+ 1 ) 

7-1 



+ const. 


(23) 


An initial condition with a shock and rarefaction emanating from the origin is constructed 
by connecting left values ( U L , V L ) to middle values (U 0 , V 0 ) with a rarefaction. The middle 
values are connected to the right values ( U T , V r ) with a shock. In our experiment the initial 
condition is given by (20), (21) where 


U L = 1.4709, U R = 2.5000, V L = 1.0388, Vr = 0.8050. 

The middle value of the solution between the shock and rarefaction is (U OJ V 0 ) = (1.973, 
1.356). We expect the viscous perturbation to have little or no effect on the speed at which 
shocks and rarefactions travel; thus, we shall compare the viscous solutions to the exact 
solution of the inviscid problem given above. 

In Fig. 6 the exact solution to the inviscid problem is compared with the computed 
solution, where Ax = .02 and e = .01. This plot shows that the shock speed and rarefaction 
propagation are nearly the same for both solutions. Greater detail is shown for the shock 
layer in Fig. 7. Here we use S x = 2e. This plot demonstrates that the method has the correct 
shock speed and that the solutions are converging to the inviscid solution as e j 0. Detail of 
the rarefaction for the same runs is shown in Fig. 8. Here the offset to the location of the 
rarefaction is Ax/2 and can be attributed to the initial condition for the parabolic layer. 

This is a domain decomposition method. The internal-layer domain is detected by com- 
paring the second derivate to tolerance of conste. The jump in internal-layer subdomain 
boundary as depicted in Fig. 9 is caused by the steady smoothing of the rarefaction until 
it is no longer detected by the residual or second derivative test. Also for t > 0.2 a simple 
shock is detected and can be solved by the ODE layer method. 
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Inviscid U Comparison With Numeric 
Figure 6: Viscous numeric compared with inviscid exact, i ■ 



Convergence to Invisicid Shock 


Figure 7: Viscous numeric compared with inviscid exact, t 




Figure 9: Internal-layer subdomain for a shock-rarefaction interaction. 
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